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ABSTRACT 



The feasibility of applying the principles of matched 
field processing to ocean acoustic tomography were studied 
under various conditions of ambient noise. Several likeli- 
hood estimators were examined (e.g., Bucker, Bartlett, 
Maximum Likelihood, etc.). Simulations were initially 
conducted for the simple case wherein only one parameter of 
the medium was unknown (e.g., SOFAR axis depth, surface 
sound speed, position of a single acoustic front) . The 
method was then applied to the more realistic problem of 
locating the boundaries of an eddy in the ocean. For 
moderate signal-to-noise ratios, all the estimators were 
shown to be able to solve the problem, albeit with different 
efficiencies. For low signal-to-noise ratios, the MLM 
scheme proved to be the most reliable especially when a 
highly correlated ambient noise was present. In all cases, 
computer simulations illustrated that mismatching may occur 
when the parameterization of the medium is poorly 
approximated. Mismatching leads to a decrease in the 
efficiency of the estimators but it may be still possible to 
correctly estimate the environmental characteristics. 
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INTRODUCTION 



I . 

As sound waves propagate through the ocean, the complex 
acoustic pressure field which is generated by the source 
depends mainly on the path followed by the acoustic rays and 
the sound speed along this particular path. Due to this 
close relationship between the sound speed field and the 
acoustic pressure field, an attempt may be made to estimate 
the range-dependent sound speed profile (SSP) between a 
fixed source and an array of receivers. The characteriza- 
tion of the SSP from acoustical measurements generally 
involves inverse techniques in order to infer the acoustic 
properties of the medium from the pressure field measured at 
the receivers. 

Due to the complexity of the ocean the inverse problem 
is most often non-linear and underdetermined. Classicial 
acoustic tomography solves this problem by linearization. 
The tomographic method is able to estimate the perturbations 
of the sound speed field by comparing the measured travel 
times of particular rays with those computed numerically 
from a reference sound speed field and a raytrace or a 
normal mode algorithm (Munk and Wunsch, 1978). The 
procedure provides maps of the perturbations in the sound 
speed field and indicates how different the actual field is 
from the one used as a reference. If a large enough number 
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of ray paths are used in the computations, the spatial 
resolution in the map of the sound speed field may be better 
than the one obtained from discrete CTD measurements (Howe, 
1986) . 

Matched field processing is a different type of inverse 
method which was first proposed as a method to locate an 
acoustic source in the ocean. The principle is to compare 
the measured complex acoustic pressures at a vertical array 
with those computed from an acoustic model using various 
positions of the target source (Bucker, 1976) . The 
procedure generates a function which is a measure of the 
likelihood between the actual acoustic pressure field 
created by the source (unknown position) and a replica 
pressure field generated from an estimate of the source 
location . 

This study is an attempt to use matched field processing 
as an alternate tool to solve the inverse problem in 
acoustic tomography. Given the position of the source and 
the receivers, matched field processing is employed to 
compare the true complex acoustic pressures at the receiving 
array with the ones computed from an acoustic propagation 
model and various sound speed fields. Computer simulations 
are used to demonstrate the performance of various 
estimators under different signal-to-noise ratios and noise 
correlation matrix structures. 
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Chapter II provides a theoretical presentation of the 
likelihood estimators which are used in this study. It also 
illustrates how the noise is modeled and added to the 
simulated data. The simulations are shown in Chapter III 
for various conditions of noise. The first simulations deal 
with the simple case where only one parameter of the medium 
is unknown (e.g., SOFAR axis depth, surface sound speed, 
single frontal boundary) and where the noise is absent. The 
next cases are applied to the localization of an eddy in a 
noisy medium. Situations of both spatially uncorrelated 
noise and correlated noise were examined. Comparison of the 
estimators is provided in Chapter IV; the spreading of each 
likelihood function about the true value is examined in more 
detail under several conditions of noise power and noise 
correlation. Also discussed is the problem of incomplete or 
poor knowledge of the other environmental parameters, e.g., 
incorrectly specifying the bottom absorption property, and 
how this may introduce inconsistency in the procedure and 
lead to a decrease in the efficiency of the likelihood 
functions . 
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II . 



ACOUSTIC TOMOGRAPHY AND MATCHED FIELDS 



A. PRINCIPLES OF MATCHED FIELD PROCESSING 

Classical beamforming for plane waves is obtained by 
measuring the maximum likelihood between the actual value of 
the complex signal at each hydrophone and the values 
computed from an expected bearing (Ziomek, 1985) . 

In a similar fashion, the distance to a target in the 
near field can be estimated by comparing the actual values 
with those computed for different distances. The range is 
assumed to be correct when both sets of values correspond to 
the same wave front curvature (Ziomek, 1985). 

Matched field processing has been used traditionally to 
find the location of an acoustic source in a shallow water 
environment. The general principle is to store the values 
of the received signal (amplitudes and phase) at each 
element of the array and then compare them with theoretical 
values computed for different possible positions of the 
target source. The true location of the emitter is 
determined when both fields match (Bucker, 1976; Baggeroer, 
Kuperman and Schmidt, 1988). 

Different criteria may be used to measure the likelihood 
or degree of matching. Each one generates a different 
function which is generally well adapted for a particular 
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type of noise. Matched field detection is consistent when 
the medium is completely determined. 

Matched field tomography deals with the inverse problem. 
Given that the source location is known, the purpose of the 
procedure is to estimate the medium characteristics, 
particularly the range-dependent sound speed profile. In 
this case, the replica or estimated fields are built from 
many sound speed profiles and one tries to match them with 
the measured one. 

B. THEORY IN NOISE-FREE CONDITIONS 
1 . Bucker Method 

According to Bucker (1976), the following "detection 
factor" may be used as a measure of the difference between 
the exact and the replica fields: 

NR NR 

/ l <KS • KR., 

BUCK = : =1 k=: ^j Z Z_ (2.1) 

F 

where terms are defined as follows: 

KR = spatial autocorrelation matrix of one replica 
field , 

KS = spatial autocorrelation matrix of the actual 
field, 

NR = number of hydrophones in the array, 

F = scaling factor to insure a result between 
0 and 1 , 
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<> = time average used to remove the component due 
to the noise, 

• = complex conjugate 

The matrices KS and KR are defined by 



KS jk “ EjEk* 



( 2 . 2 ) 



KRjk ~ E'jE'k* 



(2.3) 



where pj and p'j denote, respectively, the complex envelope 
of the acoustic pressure and the replica pressure at 
hydrophone j ; similarly for hydrophone k. As demonstrated 
by Bucker (1976), it is convenient to define the correlation 
matrices from the complex envelopes of the signal because 
the rapidly-varying time component is removed from the 
computations. These complex envelopes are easily obtained 
by processing the incoming signal through a classical 
quadrature demodulator (multiplication by a sine wave 
followed by a low-pass filter) . 

In the absence of noise, the time average is 
unnecessary (KS is time independent) and the normalized 
detection factor becomes: 



BUCK = 



NR NR 

I I ^ 

j=l k=l/j 



.. KR., 
jk* 3k 



NR 



NR 

( y 

j=l k=l/j 



KR 



... - KR., 
3 k 3 k 



1/2 ^ 

.) V2 ( I , 

j=l k=f/j 



NR 

V KS. k KS. k .) 



1/2 



(2.4) 
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This factor is similar to the classical correlation 
coefficient of two random variables. The expression given 
above does not use the diagonal elements of the matrices and 
can be interpreted as the output of a regular beamformer. 
Its value is one in the case of complete equality of the 
fields (KR = KS) . We note also that the Bucker detection 
factor is one when the matrices KR and KS are proportional, 
as a consequence of the Schwarz inequality. 

2 . Heitmever Method 

Heitmeyer (1984) defined another detection factor 
which he called the "source location ambiguity function." 
Its value is given by the following expression: 



HEIT 



NR 





NR 



NR 

V 

n=l 



En 



■ 2 



(2.5) 



and may be rewritten in terms of phases and magnitudes: 



HEIT 



NR 

I l Pi 

n=l 



■j ( ^n^n' ) 



n - 1 - n 



NR 2 
NF. I p n ' 
n=l 



( 2 . 6 ) 



It can easily be seen that the function is 
unnormalized. From the inequality: 
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HEIT < 



(2.7) 



NR 9 

I < Pn P n'» 
n-1 

NR 9 
NR l P n ' 
n=l” n 



NR NR 0 

^. p n ^ P n' 
n=l n=l 

NR 2 

NR T p ' 
t n 
n=l 



an upper limit is found: 





NR 2 

HEIT < r±r j p 
- NR S*n 
n=l 




(2.8) 


The 


ambiguity function is always 


bounded 


by a 


quantity that may be considered as the 


average 


power 


detected at 


each hydrophone. 






When 


the replica and the measured 


fields 


match 



exactly, the expression reduces to the following: 

E 1 n = En (2.9) 



which yields 



HEIT 



! NR J2 

I r ^ 

i 2 P n 

n=l n I 



NR 



NR l p, 



n=l 



n 



_1 

NR 



NR 




2 



( 2 . 10 ) 



The inequality demonstrated previously becomes an identity 
when the actual and replica fields are identical. 
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In order to obtain a normalized ambiguity function, 



we divide Equation (2.6) by this upper limit, 



HEIT = 



NR 

I I ,P n p n' exp'jtVV” 

n=l 

NR ~ NR 0 

Ip' Ip 

n=l n n=l n 



( 2 . 11 ) 



3 . Relation between Bucker and Heitmever Methods 

Using the definition of the spatial autocorrelation 
matrix, Equation (2.4) may be written: 



NR NR 



BUCK = 



: Ej- £j- ?k 'Ek 

1=1 k=l/i J 



NR NR •, n NR NR 1/2 

( T I 7 ( I L P^-Pi , 2 k ) 

j=l k=l/j ^ ^ 11 ^ j=l k=l/j ^ D ^ 



( 2 . 12 ) 



If we allow the subscripts j and k to be equal, this 
expression becomes: 



r:R mr 






NR NR 






1/2 



hence , 
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BUCK 



(2.14) 



NR 



NR 

r 



_ j 



JPj'Pj* .^ p kPk’ 

1=1 J J k=l 



NR 

l Pj 
j=l 11 



' 2 



NR 



k=l 



or, by changing the subscripts, 



NR 



BUCK = 



; n=l 



p p ' • 



NR «■) NR q 

r , Z r Z 
> p 1 / p 

S'u-i n 

n=l n=l 



(2.15) 



which is exactly the normalized ambiguity function defined 
by Heitmeyer in Equation (2.11). 

More generally, by developing the expressions of 
both functions, we can derive the following relation between 
the Bucker and Heitmeyer definitions: 



BUCK 



HEIT 



NR NR 

I l lEiEj 

i=l j=l 1 3 

NR 

2 H p. p. ' 

i=l 1 - 1 



, 1 2 



(2 . 16) 



This relation is not a simple proportionality ratio 
because it changes with the replica fields p'j. 

4 . Center of Gravity Method 

Other likelihood estimations can be developed using 
any function which has a maximum in case of perfect 
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matching. The center of gravity method is a procedure which 
solves the problem from a more mathematical viewpoint. 

In the complex plane, vectors P and P' define two 
sets of NR points, where NR represents the number of 
receivers of the array. The components of the points are 
the real and imaginary parts of the complex acoustic 
pressures : 

set P is composed of points (p-^ cos ii,Pi sin ( j_ ) 
set P' is composed of points ( p ' cos sin v '^) 



In order to compute a detection factor, we calculate 
the euclidean distance between the centers of gravity G and 
G' of both sets (Figure 2.1). This distance is inversely 
related to the likelihood of the fields. 

With a sum of weights of 1, the points G and G' are 
given by their coordinates: 



X 



g 



_1 

NR 



NR 



P 

1 



cos 



'n 



Y 



g 



NR 



NR 



o 



n=l 



n 



sin 



n 



(2.17) 



(2.18) 



and 



X' 



g 



1 NR 

NR -, P n 
n=l 



cos <j>. 



n 



Y ' 



g 



NR 



NR 



l P 

n=l 



n 



sin <t 



n 



(2.19) 



( 2 . 20 ) 
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The distance becomes: 



D = ((Xg-Xg') 2 + (Yg-Yg ') 2 ) 1/2 



( 2 . 21 ) 



or 




( 2 . 22 ) 



We then normalize the quantity in order to obtain 
unity for complete matching. 



DN = 1 - D/D max (2.23) 

where D max is the largest unnormalized distance among all 
the replicas. 

DN is only an estimation of likelihood. Although it 
is possible for two different sets of points to have the 
same center of gravity, if they are concentric, this does 
not occur in the simulations and the method keeps its 
consistency. We will see later that this distance function 
may lead to high secondary lobes and thus is not always 
reliable . 

C. THEORY IN PRESENCE OF NOISE 

Although it is possible to use the former expressions 
when noise is present which contaminates the signal, the 
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following two functions are more specifically suited for use 
in the presence of noise. 

Johnson (1982) previously demonstrated the equivalence 
between the problem of bearing determination and the 
estimation of the spectrum of a signal. Due to this 
similarity, all modern spectral estimation algorithms apply 
equally to beamforming and matched field processing. 

Although many functions may be used, as for example, 
MUSIC (Schmidt, 1981) or linear predictor (Johnson, 1982), 
we will particularly emphasize the Bartlett and Maximum 
Likelihood parameters which are two powerful estimators in 
target location problems. These estimators are especially 
useful in noisy conditions, because Baggeroer and his 
colleagues (1988) showed that they reduce to the Bucker or 
Heitmeyer structures when the noise is absent. 

1 . Bartlett Method 

This method comes directly from spectral estimation 
theory. Baggeroer, Kuperman and Schmidt (1988) demonstrated 
that the power output of a Bartlett beamformer could be 
written in the following quadratic form: 



BART = W KT W (2.24) 

where W represents the normalized velocity potential vector 
of the replica field and KT is the total spatial correlation 
matrix of the signal embedded in noise. 
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Due to the proportionality between the velocity 
potential and the pressure field, the equivalent expression 
will be used: 

BART = (P'/IP'D* K T (P'/|P'|) (2.25) 

where P' is the complex acoustic pressure vector of the 
replica at the array. Under the condition of perfect 
matching and no noise, 

NR 2 2 

BART = ) Id. I = j P ' (2 . 26) 

i=l ^ 

which is the summation of all signal powers among the 
hydrophones . 

In order to normalize the function, we will divide 
the estimator by its largest value. For comparison between 
the different methods, we will focus on the width of the 
main lobe rather than its absolute value. 

If and ni denote, respectively, the signal and 

the noise pressure at hydrophone i, the spatial correlation 
matrix has the value: 

KTij = E( (Bj+m) (Ej+nj) ’ ) (2.27) 

where E( ) denotes the expectation operation. 
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When the signal and the noise are uncorrelated, the 
matrix reduces to the simple sum of signal and noise 



matrices : 



KTij = E(piPj-) + E(niQj-) 



(2.28) 




(2.29) 



These matrices are hermitian and at least 
semidef inite . 

2 . Maximum Likelihood Method 

In spectral estimation, the Maximum Likelihood 
method, also called Capon's method or the minimum variance 
method, is used to compte the power spectral density of a 
random process (Kay, 1988). Its expression is given by: 



P MV (f) = ( e H R^" 1 e)" 1 



(2.30) 



where : 



f = frequency, 



R xx = time correlation matrix of the process, 
e = vector whose i th component is e ^ 2 f, 



H = transposition of the conjugate matrix. 



In a similar way the output of a Maximum Likelihood 



beamformer is defined: 
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(2.31) 



MLM = (W- KT -1 W) -1 

where W and KT have previously been defined. 

Following the procedure of Equation (2.25), Equation 
(2.31) is modified: 

—r -V -> -t 

MLM = ((P'VIP'I) KT -1 (P'/IP'I)) -1 (2.32) 

In the absence of noise, KT reduces to the spatial 
autocorrelation of signal only, KS (Equation (2.29)). As 
can be seen when the array is composed of two hydrophones, 
the matrix is generally singular and has no inverse. The 
calculation is made possible by adding a small amount of 
noise to the diagonal. As with the Bartlett estimator, the 
Maximum Likelihood factor will be divided by its largest 
value for normalization. 
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Figure 2.1 Distance Between Two Complex Sets of Points 
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Ill . SIMULATION OF ACOUSTIC TOMOGRAPHY 



Matched field processing has the same form as classical 
beamforming. However, instead of comparing the actual field 
vector with a plane wave replica vector, we will try to 
match the actual vector with a vector computed from an 
acoustic propagation model. Measured data will also be 
simulated with the same code, then embedded or not in noise 
depending on the scenario under investigation. 

A. PROCEDURE 

1 . Description of the Simulation 

The receiver is modeled as a vertical array and is 
assumed to be composed of 20 hydrophones evenly distributed 
between 550 m and 1500 m. The source is located 100 km from 
the receiver at a depth of 1000 m. It emits a pure sine 
wave (tonal) centered at 100 Hz. Due to the inherent 
limitation of vertical angles in the parabolic approximation 
(Ziomek, 1985) , the transmitter was selected to have a 
beamwidth of 40°. The bottom is 5000 m deep and is assumed 
to be flat and fully absorbing. This choice was made to 
speed up the calculations and is not a restrictive 
assumption. It assumes all the energy propagates by water- 
borne paths. 
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2 . Acoustic Propagation Model 

Several models could have been used to simulate the 
acoustic fields, for example, normal mode theory or the 
parabolic equation (PE) . The PE model was used because it 
is more suitable for a deep water simulation; for example, 
an ocean bottom of 5000 m allows almost 670 propagating 
modes at 100 Hz and would have been computationally 
intensive using normal mode theory. 

The classical PE approximation with split-step 
Fourier transform (Coppens, 1982) was available in the 
Environmental Acoustic Research Group package of models 
resident at NPS. The source code was slightly modified to 
save the complex acoustic pressures at the hydrophones in a 
file. The measured and the replica complex pressure fields 
were then stored in order to run the simulation programs. 

3 . Simulation of Noise 

Noise was added to the measured data in order to 
produce a realistic problem and to study the behavior of the 
estimators in different environments. Following the 
procedure described by Porter, Dicus and Fizell (1987), 
noise was introduced by the mean of its spatial correlation 
matrix. This procedure is better than just altering the 
data with random noise because the probability density 
function of the noise is difficult to estimate. Moreover, 
all the estimators considered were written in terms of 
correlation matrices. 
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Ambient noise falls in two categories: 

- uncorrelated noise. Its correlation matrix is 

proportional to the identity matrix and the 
proportionality factor is an indicator of the noise 
power. 

- correlated noise. In this case the matrix has non zero 
terms outside the diagonal, but is nevertheless an 
hermitian matrix. 

Several attempts to measure the coherence of ambient 
noise in the ocean have been conducted during the past 
years. One of them was made from the Trident Vertical Array 
and is described by Urick (1984) . Figure 3.1 depicts the 
results of this study and has been used to generate a model 
of the noise correlation matrix. 

The matrix was modeled in the following way: 

KNij = a 2 e - I i- 3 I (3.1) 

where j 2 depends on the noise power and 3 is a factor which 
indicates how fast the coherence falls off outside the 
diagonal. The larger 8 is, the more uncorrelated the noise 
is, i.e., its spatial correlation scale becomes shorter. 

Due to the spacing between the receivers in the 
array and the frequency (100 Hz) used in the simulation, 
Figure 3.1 shows that 6= 1.7 is consistent with the 
observed ambient noise correlation. 

Although the spatial correlation matrix of the 
noise, KN, is generally a complex hermitian matrix, the 
analytic modeling shown in Equation (3.1) describes a real 
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symmetric matrix. As explained by Cox (1973), this 
approximation is valid in the special case of zero time 
delay, i.e., when it is assumed that the noise is in phase 
among all the hydrophones of the vertical array. This 
assumption is relatively consistent for low frequency noise. 
Below 150 Hz the noise is principally due to distant 
shipping and arrives mainly from the horizontal. 

B. ONE DIMENSIONAL PROBLEM (NOISE-FREE CONDITIONS) 

In the following simulation, the shape of the sound 
speed profile is the only unknown. If the profile is 
digitized in 1 m intervals, then for a water depth of 5000 
m, one would have to determine 5000 values to match the 
complete sound speed profile. Such a procedure would lead 
to an unmanageable number of computations, especially if one 
tries to match a large number of replica fields. However, 
as we are only interested in demonstrating the feasibility 
of the procedure, we will begin by studying the simple cases 
where only one or two points of the sound speed profile are 
unknown . 

1 . Determination of a SOFAR Axis Depth 

We initially start with a bi-gradient sound speed 
profile having a sound speed minimum at 1000 m (Figure 3.2). 
Two replica profiles are considered wherein the SOFAR axis 
depth is altered by +/- 200 m (step 10 m) . Note that 
because the surface and bottom sound speeds are unchanged. 
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the change in axial depth results in a change in the 
gradients of both the upper and lower segments of the SSP. 

Three estimation techniques, the Bucker, Heitmeyer, 
and center of gravity methods are utilized to determine the 
true depth of the sound speed minimum. These estimators are 
computed from Equations (2.4), (2.11) and (2.23). The 
estimated depth of the SOFAR axis is found when an estimator 
shows a peak with a detection factor of 1. The results are 
shown in Figures 3.3, 3.4 and 3.5; each estimator 
demonstrates a different behavior. 

The Bucker detection factor and the Heitmeyer 
ambiguity function both indicate a maximum at the true 
location of 1000 m. However a strong side lobe, centered at 
1040 m, indicates these two detection factors are not robust 
enough to provide an unambiguous selection of the SOFAR axis 
depth. In addition, strong secondary side lobes are also 
present . 

The center of gravity method (Figure 3.5) proves to 
be a better estimator for this situation. The main lobe is 
much narrower and no other lobes exist. For a noise-free 
ocean, this is the best estimator among the three to solve 
this particular problem. 

2 . Determination of Surface Sound Speed 

In order to mimic typical seasonal or spatial 
changes in the SSP, the surface sound speed was permitted to 
fluctuate by +/- 5 m/s about a mean value (step 1 m/s). For 
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this situation the SOFAR axial depth was fixed at 1000 m. 
Hence, only the upper gradient changes as seen in Figure 
3.6. We seek an estimator that will match the true surface 
sound speed (solid line in Figure 3.6). Using the same 
three estimation techniques as above, the results are shown 
in Figures 3.7, 3.8 and 3.9. For this situation, we observe 
a nearly identical behavior of the Bucker and the Heitmeyer 
functions, both of which have moderate side lobes at about 
0.7. As before, the center of gravity estimator remains the 
best without any ambiguity due to the presence of secondary 
lobes . 



3 . Determination of an Acoustic Frontal Boundary 



To 


model 


the 


presence 


of an acoustic 


front two 


different 


sound 


speed 


profiles 


are introduced, 


one 50 km 



from the other. The parabolic equation model was utilized 
in this range-dependent problem with the position of the 
front (i.e., the range at which the second SSP is 
encountered) allowed to vary by +/- 2 0 km about the true 
position (step 1 km). Figure 3.10 provides an illustration 
of the SSP setup. The Bucker and the Heitmeyer methods 

correctly solve this problem with relatively narrow main 
lobes, as shown in Figures 3.11 and 3.12. However, the 
center of gravity method does not perform as well (Figure 
3.13). Although it is able to locate the correct value, 
significant sidelobes are present which could lead to an 
ambiguity if too small a detection threshold were chosen. 
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This last procedure is obviously not suitable for these 
conditions. 

4 . Comments on the Estimators 

In order to completely appreciate the consistency of 
each of the previous estimators for the case where only one 
factor is unknown, it is important to study their response 
to a variety of configurations. 

a. Influence of the Number of Hydrophones 

In the problem of determining the depth of the 
SOFAR axis, the array was composed of 20 hydrophones. It is 
possible to run the same simulation by using only a fraction 
of the receivers. Plots of the Bucker estimator for four 
different numbers of hydrophones are shown in Figure 3.14. 
Although the difference is small when the number of 
hydrophones is reduced from 2 0 to five, the output of an 
array composed only of two receivers changes drastically. 
When the number of hydrophones is this small, the side lobes 
may have an amplitude of the same order as the main lobe, a 
situation which leads to full ambiguity. The number of 
hydrophones is thus an important parameter which must always 
be more than some minimum value. This value is 

unfortunately dependent upon the actual problem and the 
depth of the array relative to the axial depth. Receivers 
which do not intercept much of the acoustic energy can be 
easily omitted but those which contain significant amplitude 
and phase information should be retained. 
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b. Influence of the Frequency 

Typically as a result of beamf orming , the main 
lobe becomes narrower as the frequency of the array 
increases. The same phenomenon can be observed in Figure 
3.15, wherein the width of the estimator peak is also a 
function of the frequency. However, a trade-off exists 
between the desired resolution (width of the detection peak) 
and the computation time of the PE model which increases 
rapidly with frequency. Also if higher frequencies (kilo- 
Hertz range) were used, the signal would be limited by 
absorption which would result in lower signal-to-noise 
ratios . 

c. Importance of Array Position 

The depth of the array is of minimal importance 
when working in shallow water because the entire water 
column is nearly insonified. Such is not the case in deep 
water where shadow zones exist with relatively low signal 
levels. Based on the depth of the source, a first guess of 
the depth to position the array would be to place it where 
the signal may be expected to occur with a high level, for 
example, in the vicinity of the SOFAR axis or near a 
convergence zone. The choice will obviously be dependent on 
the profile shape. 
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C. TWO PARAMETER PROBLEM (NOISY CONDITIONS) 

1 . Localization of an Eddv 

The previous section has shown that the estimators 
are generally able to find the correct value in the case of 
a simple unknown and a noise-free medium. The same kind of 
simulation may be run when two parameters are to be 
determined. To test the ability of the various estimators 
to deal with a two dimensional problem, we will examine 
their ability to locate an eddy assumed to be present 
between a source and a receiving array. The sound speed 
profiles inside and outside the eddy are known. Thus the 
only unknowns are the borders of this perturbation of the 
sound field. An eddy 20 km in diameter is positioned 40 km 
to 60 km from the source. The replica fields are computed 
by scanning the limits from 35 km and 45 km for the border 
closest to the source, and from 55 km to 65 km for the 
farther boundary. Replications at 1 km interval were made. 
We will thus try to match the simulated measured data with 
121 replica fields. Figure 3.16 presents the true location 
of the eddy in this simulation. 

Figure 3.17 shows the disposition of source and 
array with regard to the energy field for the true location 
of the eddy. The array lies on the SOFAR axis, almost 30 km 
beyond a convergence zone. From this plot, we can expect a 
high level of signal from the channel propagation and large 
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differences in phase due to multipath propagation (RR 
acoustic rays) . 

Before processing, the measured data are imbedded in 
noise. This noise is introduced by the spatial autocorrela- 
tion matrix described earlier. The principal assumption of 
this simulation is that all correlation matrices are 
completely known. For an actual situation, this may not be 
true but it is still possible to estimate the total matrix 
of noise from the set of measured data. 

Because of its close similarity to the Bucker 
detection factor, the Heitmeyer function will be omitted 
from further analysis. The Bartlett and the Maximum 
Likelihood functions are introduced for these simulations, 
and it will later be seen that these two estimators are well 
suited for conditions where the signal-to-noise ratio is 
low . 

2 . Signal-to-Noise Ratio 

The signal-to-noise ratio is a parameter which 
depends on the relative powers of the signal and the noise 
at the array. Following the procedure of Ziomek (1985), the 
signal-to-noise ratio can be written in terms of the noise- 
free signal and the spatial autocorrelation matrix of the 
noise : 



SNR = 20 Log 10 



NR NR 

& 

NR NR 

y kn. . 

i-1 jh « 



(3.2) 
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The numerator of this expression may be interpreted 
as the summation of all the elements of the signal spatial 
correlation matrix. Similarly, the denominator represents 
the sum of the entries of the noise correlation matrix. 

The SNR was computed for several values of the noise 
power, a 2 , and the correlation parameter, S, that generate 
new values of KN-^j (Equation (3.1)). The results are 
presented in Table 3.1. The expression above shows that the 
SNR is reduced when the denominator of the argument 











NR 


NR 




increases , 


i . e . , 


when the 


double 


summation 

i=l 


KNij 

3=1 


is 


large due 


to a 


significant 


increase in the power of 


the 


noise, ~ 2 , 


or a 


slow decay 


in the 


correlation 


between 


the 



hydrophones, . By using the results presented in Figure 
3.1 and the analytic expression of the noise matrix given in 
Equation (3.1), the correlation of the noise will be 
considered high when the factor 6 is less than 0.57 and low 
when it exceeds 2.0. 



TABLE 3 . 1 

SIGNAL TO NOISE RATIO AS A FUNCTION OF ~ 2 AND 2 
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In the following sections, the performance of 
various estimators will be examined for several conditions, 
including different cases of noise power and noise 
correlation. 

3 . Bucker Detection Factor 

A simulation using the Bucker detection factor was 
run for the case of eddy localization under both noise-free 
and noisy conditions. 

a. Noise-free Conditions 

The simulation was run by setting the power of 
the noise, to zero. The autocorrelation matrix of the 
noise is then just the null matrix and the total matrix 
reduces to one of signal only, as indicated by Equations 
(3.1) and (2.29). A three dimensional plot and a contour 
plot of the detection factor are shown in Figure 3.18. The 
estimator is represented by a surface which has only a 
single maximum positioned at the correct location of the 
eddy. This suggests that the Bucker method is able to 
determine the true location of the eddy in a noise-free 
environment . 

An interesting feature of the plot is the 
symmetry that exists around both diagonals of the contour. 
Moving on the principal diagonal, along the line: 



Y = X + 20 



(3.3) 
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is equivalent to displacing an eddy with a constant diameter 



of 20 km. The small intervals between the contour lines 
along this path indicate that the location of the eddy can 
be determined with reasonable accuracy, once we know its 
diameter . 



large value of 3 . In this case, the correlation falls off 
rapidly on either side of the noise matrix diagonal. For a 
large enough 3 , the noise field at one hydrophone is 
completely dissimilar to that at another hydrophone and the 
noise matrix becomes diagonal. Simulations were run with 
t= 10 and different values of o^. All yielded the same 
results as in Figure 3.19, which is seen to be identical to 
Figure 3.18. This similarity may be explained by recalling 
the definition of the Bucker detection factor when noise is 
present : 



b. Case of Uncorrelated Noise 



Uncorrelated noise is generated by choosing a 



MR NR 



BUCK 




(2.4) 



where the total correlation matrix is given by 



KTj k = KS jk + KN jk 



(2 . 29) 
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or, 



KT jk = KS jk + - 2 exp(-i|j-k|) (3.4) 

As can be seen, this expression only uses the 
cross terms of the matrix. When 8 is large enough, the 
second term in the right hand side of Equation (2.29) is 
almost negligible for every value of the noise power, " 2 ; 
thus the cross terms of the total matrix KT reduce to the 
cross terms of the correlation matrix of the signal KS . 

KT j k ' KS j k , j 4 k (3.5) 

By changing the value of - 2 , the power of the 
noise is modified, but the new diagonal terms do not play a 
role in the calculations. The Bucker detection factor is 
thus insensitive to perfectly uncorrelated noise; in this 
case, the performance is exactly identical to the one in 
noise-free conditions. 

c. Case of Uncorrelated Noise 

Any combination of noise power, o^ f and noise 
correlation, g, yields a different pattern of the detection 
factor. In cases for which the spatial correlation of the 
noise is high, the Bucker method may still be able to 
maximize the detection factor at the correct location, but 
the absolute value of the peak will decrease as the 
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ambiguity surface becomes flatter. Figure 3.20 illustrates 
this type of behavior. 

Whenever both 3 and o 2 generate a low signal-to- 
noise ratio (large power, o 2 , or small correlation 
parameter, 8 ) , the procedure fails and the localization of 
the eddy becomes impossible (see Figure 3.21). We will 
quantify the effects of a 2 and g on localization below. 

4 . Bartlett Estimator 

The same simulations as above were run using the 
Bartlett estimator under noise-free and noisy conditions. 

a. Noise-free Medium 

In a generic noise-free environment, the 
performance of the Bartlett estimator is similar to the 
Bucker detection factor, as shown in Figure 3.22. As 
expected, the same symmetry along the diagonals is still 
present . 

b. Correlated Noise 

The performance of the Bartlett estimator 
changes significantly as ^ 2 and vary. Figure 3.23 

provides an example of the plot for c 2 = 10 -1 ^ and 8 = 1.7, 
where the true location is found. Figures 3.24 and 3.25 
illustrate failures of the method due to a weak signal-to- 
noise ratio brought about by strong noise and highly 
correlated noise, respectively. It will later be 

established that the usual characteristics of the noise in a 
deep ocean do not generally lead to this kind of ambiguity. 
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5 . Maximum Likelihood Method 

The MLM estimator was calculated using the same 
noise conditions as for the Bartlett function. An 
examination of noise-free conditions is not possible because 
of the singularity of the correlation matrix. Equation 
(2.31) shows that the expression of the Maximum Likelihood 
estimator requires the calculation of the inverse matrix 
KT -1 . 



MLM = (W* KT -1 W) -1 (2.31) 

When noise is absent, the matrix KT reduces to the 
correlation matrix of the signal KS which is generally 
singular. 

a. Slightly Correlated Noise 

When 2 = 10 -16 and = 1.7, the method gives 
better results than the Bartlett estimator with relatively 
low side lobes (compare Figure 3.26 with Figure 3.23). 

b. Strongly Correlated Noise 

With a more highly correlated ambient noise 
(s = 0.17), as depicted in Figure 3.27, it is still possible 
to obtain a correct location of the eddy. By comparing this 
plot with Figure 3.25, we note that the Bartlett estimator 
was unsuccessful in this case. Nevertheless, even for the 
MLM technique, a very low signal-to-noise ratio will result 
in a failure as shown in Figure 3.28. 
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Figure 3.1 Spatial Correlation of Ambient Noise for 

Frequencies 50-100 Hz ( ) and 200-400 Hz 
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Figure 3.2 Actual and Extreme Replicas of the SSP used 
in the Simulation to Determine SOFAR Axis 
Depth 
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Figure 3 . 3 



SOFAR Axis Depth Determination, 
Detection Factor 
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Figure 3 . 4 



SOFAR Axis Depth Determination. 
Ambiguity Function 
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Figure 3.5 SOFAR Axis Depth Determination, 
of Gravity Method 
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Figure 3.6 Actual and Extreme Replicas of the SSP used 
in the Simulation to Determine the Surface 
Sound Speed 
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Figure 3.7 Surface Sound Speed Determination 
Detection Factor 
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Figure 3.8 Surface Sound Speed Determination 
Heitmeyer Ambiguity Function 
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Figure 3.9 Surface Sound Speed Determination, 
of Gravity Method 
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Figure 3.10 Simulation to Determine the True Location 
of an Acoustic Front 
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Figure 3.11 Acoustic Front Determination. Bucker 
- Detection Factor 
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Figure 3.12 Acoustic Front Determination. Heitmeyer 
Ambiguity Function 
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Figure 3.13 



Acoustic Front Determination, 
of Gravity Method 
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SOFAR Axis Depth Determination for Various 
Numbers of Hydrophones in the Array 
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Figure 3.14 
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Figure . i- Surface Sound Speed Determination. 

Influence of the Frequency of the Source 
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Figure 3.16 Eddy Localization. True Position ( 
of the Eddy Boundaries 
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Figure 3 . 18 



Eddy Localization. Bucker Detection 
Factor. Noise-Free Conditions 
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Figure 3.19 Eddy Localization. Bucker Detection Factor 
Illustrating Condition of Uncorrelated 
Noise and Moderate Noise Power, = lO - ^-®, 
6 = 10 
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Figure 3.20 Eddy Localization. Bucker Detection Factor 
Illustrating Condition of Highly Correlated 
Noise and Moderate Noise Power, " 2 = 10~ 17 , 
P = 0.15 
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Figure 3.21 Eddy Localization. Bucker Detection Factor 
Illustrating Condition of Highly Correlated 
Noise and Strong Noise Power, = 10 ~ 16 , 

0 = 0.15 
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Figure 3.22 Eddy Localization. Bartlett Estimator. 
Noise Free Conditions 
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Figure 3.23 Eddy Localization. Bartlett Estimator 

Illustrating Condition of Moderate Noise 
Power and Moderate Noise Correlation, 
o 2 = 10 -16 , 6 = 1.7 
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Figure 3.24 Eddy Localization. Bartlett Estimator 
Illustrating Condition of Strong Noise 
Power and Moderate Noise Correlation, 

o 2 = 10 -15 , 6 = 1.7 



57 




ctosr: limit or loor ikm) 



Figure 3.25 Eddy Localization. Bartlett Estimator 

Illustrating Condition of Moderate Noise 
Power and Strong Noise Correlation, 
o 2 = lo -16 , & = 0.17 
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Figure 3.26 Eddy Localization. MLM Estimator Illus- 
- trating Condition of Moderate Hoise Power 
and Moderate Noise Correlation, o 2 = io -16 

6 = 1.7 
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Figure 3.27 Eddy Localization. MLM Estimator Illus- 
trating Condition of Moderate Noise 
Power and Strong Noise Correlation. 

„2 = 10 - 16 , p , 0 . 17 
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Figure 3.28 Eddy Localization. MLM Estimator Illus- 
trating Condition of Strong Noise Power 
and Moderate Noise Correlation, o2 _ ^q-15^ 
8 =1.7 
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IV. ANALYSIS OF THE SIMULATION 



From the previous simulations, one sees that the 
performance of each of the estimators varies significantly 
under varying signal-to-noise ratios or source/receiver 
geometries or sound speed variations. For example, the 
center of gravity method was shown to be the best in 
locating the SOFAR axis and determining the surface sound 
speed. In contrast, this procedure was the least successful 
in the acoustic front localization problem. Therefore the 
efficiency of an estimator does not depend only on the type 
of ambient noise but also on the particular problem being 
solved. In order to continue focusing on a realistic 

problem, the eddy localization problem will be studied in 
greater detail. The condition of mismatching will be 
treated separately. 

A. COMPARATIVE STUDY OF THE ESTIMATORS 
1 . Criterion 

To facilitate comparison of the performance of the 
different methods, we will calculate the joint central 
moments of the different estimators. The joint central 
moment provides information on the spread of the detection 
factor about the mean. The smaller the moment, the better 
the estimate of the parameter of the ocean. 
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Using the definition of the central moment of a 
multiple random variable as defined by Peebles (1987): 



^nk = / /(x-X) n (y-Y) k f xy (x, y) dx dy (4.1) 



where : 



X, Y = means of the random variables X and Y, 

f X y = joint probability density function of X and Y, 

n, k = orders of the central moment. 

The spreading factor for the MLM method is defined 
as the second order central moment of the function MLM(x,y): 

S = ' MLM (x, y ) ( x-4 0 ) 2 (y-60) 2 dx dy (4.2) 

X V 

where x and y are the boundaries of the eddy we are looking 
for . 

Since we are dealing with a discrete search among 
parameter values, this pseudo variance has the form: 

NR NR 

S = l l MLM-h ( 44+i — 40) 2 ( 54 + j -60) 2 (4.3) 

i=l j=l 

where NF denotes the number of possible values of i and j. 

Equation (4.3) indicates how the estimator spreads 
around the true frontal boundary values x = 40 and y = 60. 
It is a global measure of the estimator performance, not 
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just a measure of the main lobe width. A large value is 
associated with significant spreading which indicates a poor 
performance, even if the true position is actually found. 
The value of the spreading factor will be also large when 
the ambiguity surface has a large mean value and a small 
amplitude (case of significant noise). Analogous spreading 
factors may be defined for the Bucker and the Bartlett 
functions . 

2 . Efficiency of the Estimators 

The spreading factor defined above has been computed 
for the different combinations of and ' shown in Table 
3.1. The results are presented in Table 4.1 for the Bucker, 
Bartlett and MLM methods. From this table it is possible to 
compare the efficiency of each technique in a variety of 
environments. The values of the correlation factor i have 
been chosen to stay consistent with deep water measures. 
The range of noise powers, produced signal-to-noise 
ratios between -67 dB and +21 dB. 

In order to be consistent in comparing the spreading 
factor S of the three schemes, it is convenient to modify 
the expression for the Bucker detection factor as defined by 
Equation (2.4) by dividing it by its largest value among all 
the replicas. The maximum of this new function will always 
be one in all situations of noise, as the Bartlett and the 
MLM estimators. 
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TABLE 4 . 1 



SPREADING FACTOR S OF THE BUCKER (NORMALIZED) , BARTLETT 
AND MLM ESTIMATORS AS A FUNCTION OF c 2 AND 



2 

o 



10-19 



5 *10 -19 



10" 18 



5 *10 -18 



10 - 1 7 



10 -16 



10-15 



0 . 57 

SNR=+13dB 
19 , 000 
25,924 
114 

SNR=-ldB 

19,280 

26,471 

567 

SNR=-7dB 

19,632 

27,149 

1,131 

SNR=-2 ldB 
22 ,400 
32 , 328 
5 , 499 

SNR=-27dB 
25,773 
38 , 256 
10,630 

SNR=-4 7dB 



66,998 

SNR=-67dB 



1 .00 

SNR=+16dB 
18,973 
25 , 902 
139 

SNR=+3dB 

19 . 147 

26.364 
654 

SNR=-3dB 

19 .364 
26,935 

1,382 

SNR=-17dB 
21,081 
31 , 320 
6,657 

SNR=-23dB 
23 , 188 
36,371 
12 , 732 

SNR=-43dB 
54 , 747 
85,126 

72 . 147 

SNR=-63dB 

oc 

00 

00 



1.70 

SNR=+2 OdB 

18 . 950 
25,884 

142 

SNR=+6dB 

19,031 

26,272 

708 

SNR=0dB 

19,133 

26,754 

1,411 

SNR=-14dB 

19,946 

30,463 

6,766 

SNR=-2 OdB 

20.950 
34 , 771 
12 , 883 

SNR=-39dB 

37,381 

78,315 

70,060 

SNR=-60dB 

00 

oc 

00 



2 . 00 

SNR=+2 ldB 
18 ,944 
25,880 
140 

SNR=+7dB 

19,004 

26,250 

698 

SNR=+ldB 

19,078 

26,710 

1,391 

SNR=-13DB 

19,671 

30,258 

6,670 

SNR=-19dB 
20,407 
34 , 386 
12 , 687 

SNR=-39dB 

32,723 

76,655 

69,050 

SNR=-59dB 



2 .20 

SNR=+2 ldB 
18,941 
25,877 
139 

SNR=+7dB 
19 , 090 
26,239 
691 

SNR=+ldB 
19,050 
26,688 
1, 377 

SNR=-12dB 
19 , 532 
30,153 
6,603 

SNR=-19dB 
20,131 
34 , 190 
12 , 570 

SNR=-3 8dB 
30,280 
75,805 
68 ,404 

SNR=-58dB 

X 

124,781 
124 , 584 



Note 1: The signal-to-noise ratio is indicated for each 

entry of the table, followed by a triplet of numbers which 
represent, respectively, the Bucker (normalized) , Bartlett 
and MLM spreading factors. 

Note 2: The value °° means that the estimator is unable to 

detect the true location of the eddy. 
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Table 4.1 provides a good illustration of the 
performance of the estimators under several conditions of 
noise power and correlation. Figures 4.1, 4.2 and 4.3 show 
logarithmic plots of the spreading factor S versus the noise 
power 0 2 ; each curve represents a value of the correlation 
parameter 3 . It is thus possible to determine the value of 
S for any pair of a 2 and S. As is obvious from Table 4.1, 
the comparative performance of each method is mostly a 
function of the signal-to-noise ratio of the measure. 

For SNR less than -50 dB all methods fail to detect 
the true location of the eddy. Nevertheless, we observe the 
case a 2 = 10“-*-^ and 3 = 2.2, where the Bartlett and MLM 
estimators indicate two different maxima at one. Even in 
this case, the amplitude of the functions is so small that 
the spreading factor S is very large. 

When the SNR is about -4 0 dB, the Bucker method is 
the most efficient estimator, albeit a weak one, as the 
spreading remains significant. The superiority of the 
scheme increases moreover when the correlation of the noise 
decreases . 

For other SNR and correlation values, the MLM method 
is generally the most efficient. When the SNR = 0 dB or 
greater, the advantage of the MLM scheme is obvious with 
respect to the other two functions. One also notes that the 
MLM estimator is well adapted to resolving highly correlated 
noise situations; in such cases, simulations show that the 



66 



width of the main lobe becomes narrower but that the mean 



component of the surface increases. 

B. MISMATCHING CASE 

Mismatching occurs when a parameter used in simulation 
of the replica fields has been incorrectly estimated and is 
different from the one that created the true data. In the 
eddy localization problem, this defect may be introduced in 
several ways: 

- a wrong measure of the source frequency, 

- inaccurate estimation of the source or array position, 

- inaccurate estimation of the source beamwidth, 

- insufficient knowledge of the bottom loss 
characterization, 

- oversimplification of the SSP. 

This list is not exhaustive and one must keep in mind 
that perfect matching almost never exists due to the 
impossibility of any acoustic model to solve the true 
acoustic wave equation in the real ocean. Because 

mismatching prevents a close likelihood between the actual 
and replica fields. It can be thought of as an additional 
noise which the correlation matrix has not taken into 
account. Mismatching thus results in a degraded estimation 
of the total spatial correlation matrix KT. The next 
section studies in greater detail two cases where 
mismatching is created by a change in the bottom loss 
parameterization and where the borders of the eddy are 
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smoother than what we have used previously to simulate the 
real measured data. 

1 . Change in Bottom Loss Parameterization 

All simulations were conducted with the assumption 
that the bottom was fully absorbing. We will now consider 
the bottom to be a perfectly rigid surface with total 
reflection and examine how this new treatment modifies the 
likelihood functions. 

In the absence of noise the difference in phase at 
each receiver of the array for both the perfectly reflecting 
and fully absorbing bottom conditions is depicted in Figure 
4.4. When the bottom is treated as a perfect reflector, 
perfect matching will not occur because all the replica 
fields are constructed based on the full absorption 
assumption . 

Figures 4.5 and 4.6 show the result of a simulation 
using the MLM estimator with characteristics c 2 = 10 -17 and 
6 = 1.7 to represent the situations of no mismatching and 
mismatching, respectively. When mismatching occurs, the 
amplitude of the peak decreases due to an increase in the 
mean value of the likelihood function. The secondary lobes 
also become larger. It is important to note that the 
degradation observed in Figure 4.6 does not imply that 
treating the bottom as a perfect reflector is less correct; 
it only means that the matching was done improperly; one 
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must therefore be consistent in modeling the environmental 
input variables. 

2 . Oversimplification of the Actual Eddy 

The previously described eddy localization 
simulations were run with an excessive simplification of the 
true medium. The simulations assumed that there were no 
horizontal gradients of sound speed in and outside the eddy 
boundaries, i.e., the change of SSP occurred almost 
instantly at the borders of the eddy. In an attempt to be 
more realistic, the next simulation was conducted after 
adding four intermediate SSPs between the two previously 
utilized profiles (Figure 4.7). The first intermediate SSP 
was introduced 4 km before the border of the eddy and the 
next ones added every kilometer thereafter. Actual signal 
values were generated with this smoothed barocl inicity , and 
replica fields were computed as before with the simplistic 
three-profile model. Noise was omitted from this simulation 
in order to better appreciate the effect of this 
mismatching. Results for the Bucker method are presented in 
Figure 4.8. Comparing this plot with Figure 3.18 we see 
that the main peak decreases but localization of the eddy 
remains possible. The implication of this simulation is 
that identification (location) of strong frontal boundaries, 
such as the north wall of the Gulf Stream or ice edge 
fronts, could be fairly exact but weaker, open ocean 
mesoscale eddies may pose more of a difficult problem 
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(assuming the same number of profiles are used to estimate 
the replica fields) . 

The above examples suggest that in cases where the 
signal-to-noise ratio is quite low, it is possible for 
mismatching to hide the true location of the maximum and so 
possibly lead to a failure of the procedure. As often as 
possible, mismatching must be avoided by a comprehensive 
knowledge of the parameters used in the replica fields 
calculations . 

C. COMMENTS ON THE PROCEDURE 

As we were only interested in using matched field 
processing in acoustic tomography, many simplifications have 
been made to run the simulations. Although the procedure 
seems to be applicable and efficient in most cases, it is 
necessary to test its applicability to more complex 
problems . 

1 . Environment 

It was assumed in all the simulations that only one 
or two parameters were unknown, for example, the surface 
sound speed or the eddy location. For actual oceanic 

situations, many more properties can be expected to vary in 
space and time. Extending this study to a heterogeneous 
medium is theoretically feasible but vastly increases the 
number of unknowns. 

Modeling a shallow water environment has not been 
considered as a possible mechanism to speed up the 
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calculations even though matched field processing remains 
possible in this kind of environment. Several limitations 
are apparent. A better knowledge of the bottom structure is 
required. The ambient noise is moreover quite complex in 
coastal waters and its spatial correlation matrix would be 
difficult to model. As bottom interaction is important for 
estimator calculations (see the mismatching case) , it is 
possible to consider the bottom loss as an unknown parameter 
and attempt to determine it through matched field 
processing. To examine this special case, the replica 
fields are generated using nine different bottom loss curves 
and one attempts to find the actual bottom loss curve. 
Figure 4.9 illustrates the different bottom loss curves that 
have been used to create the replica fields in shallow water 
(300 m) . The maximum of the likelihood estimator occurs 
when the replica bottom corresponds to the actual loss. 
However, since the various curves are so similar in shape 
and because the bottom loss has only a weak to moderate 
effect on the transmission loss, the correct bottom loss 
curve is not sharply defined. This implies that a correct 
specification of bottom loss for a low loss bottom is not 
required; not so for a high loss bottom. 

2 . Model Consideration 
a. Resolution 

In the typical target location problem, the 
medium parameters are generally considered constant. It is 
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only necessary to run the PE model once to compute the 
pressure field at different distances. However, when the 
SSP becomes the unknown in the tomography problem, the PE 
model must be run for each replica. If we wish to find the 
correct shape of the SSP from 0 to 300 m in the seasonal 
thermocline with a resolution of 1 m, the model needs to be 
run 10 3 00 times if the sound speed at each depth can have 
ten different values. This simple example illustrates the 



trade-off between 


resolution 


and computer 


time . 


In 


the 


problem of 


eddy 


localization 


, where the 


limits 


of 


the 


perturbation 


were 


allowed to 


vary over 11 


values , 


the 


PE 



model was run 121 times. For a determination of more than 
two unknowns, the basic theory of the estimators is still 
valid but the representation of the ambiguity surfaces 
becomes unachievable due to limitations in computer run 
time . 

b. Noise Approximation 

The correlation matrix of the ambient noise was 
modeled by a symmetric matrix decaying exponentially around 
the diagonal. This approximation is rather poor, even in 
deep water, because the power of the noise is assumed to be 
the same at each hydrophone. A better simulation would have 
been to consider the matrix of a noise which is, like the 
signal, a solution of the wave equation. Our matrix is 
symmetric even though the actual matrix of the complex noise 
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needs to be hermitian, because the cross-correlations are 
complex . 

3 . Correlation Matrix 

If the spatial correlation matrix of the noise were 
known, it could be introduced in the replica fields and we 
could deal with it as with a noise-free problem. 

In practice, it is not possible to separately 
compute the noise and the signal correlation matrices. The 
total matrix needs to be estimated from the noisy signal at 
the hydrophones. Several techniques are available, as for 
example the Fourier method recommended by Johnson (1982). 
As the matrix becomes only an estimate; we should expect a 
slight mismatching and then decreasing efficiency of the 
estimators . 
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LOG 10 (S) FOR SUCKER METHOD 




Figure 4.1 Eddy Localization. Normalized Bucker 
- Detection Factor. Spreading Factor S 
Function of o 2 for Various Conditions 
Correlation, B 
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as a 
of 



LOG 1 0 (S) FOR 8HRTLETT METHOD 




LOG 10 (51GM02 ) 



Figure 4.2 Eddy Localization. Bartlett Estimator 
Spreading Factor S as a Function of o 2 
Various Conditions of Correlation, 8 



for 
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LOG 10 (S) FOR iILM nETHOQ 




LOG 1 0 (SIGfin2) 



Figure 4.3 



Eddy Localization, MLM Estimator. 
Spreading Factor S as a Function of n 2 
for Various Conditions of Correlation, 6 
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Figure 4.4 Phase of Complex Acoustic Pressure at the 

Hydrophones for Fully Absorbing ( ) and 

Perfectly Reflecting ( ) Bottom 
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Figure 4 . 5 



Eddy Localization. MLM Estimator, 

= 10 - -*- 7 , p = 1.7. No Mismatching 
(Bottom Treated as Fully Absorbing) 
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Figure 4.6 Eddy Localization. MLM Estimator, 
o 2 = 10- 17 , 3 = 1.7. Mismatching 
(Bottom Treated as Perfectly Reflecting) 
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Figure 4.7 Intermediate SSPs ( ) used for Transition 

from Open Ocean SSP (cl(z)) to the Eddy SSP 
(c2 ( z ) ) 
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Figure 4 . 8 



Eddy Localization. Bucker Detection Factor. 
Mismatching Due to an Excessive Simplifica- 
tion of Horizontal Sound Speed Gradients 
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(b) 

Figure 4.9 Bottom Loss Determination. Bucker Detection 
Factor. Replica Bottom Loss Curves (a). 
Estimation of Bottom Loss Curve (b) 
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V. CONCLUSION 



Matched field processing has been shown to be an 
efficient way to solve the inverse problem in ocean acoustic 
tomography when the ocean can be characterized by a few 
parameters. The estimators which have been used (Bucker, 
Bartlett, Maximum Likelihood, etc.) were generally robust 
enough to find the actual sound speed field of the ocean 
under usual conditions of noise power and noise correlation. 

A. SUMMARY OF THE RESULTS 

For noise-free conditions or high signal-to-noise ratios 
(SNR) , all the estimators are able to correctly determine 
the actual unknown parameter of the medium. The Bucker 
detection factor was shown to be the best function when the 
SNR was moderate. 

The efficiency of the various estimators, illustrated by 
their spreading about the true value, decreased in cases of 
low SNR introduced by a large power or a high spatial 
correlation of the noise. The Maximum Likelihood method was 
shown to be the best scheme when the ambient noise was 
highly correlated because its spreading was less sensitive 
to the degree of spatial correlation than the other 
estimators . 

The effect of mismatching, when introduced in the 
simulations, generated a decrease in the efficiency of the 
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methods. Analysis of several degrees of mismatching 
indicated that unambiguous results can be expected from 
matched field processing provided that the parameterization 
of the medium is exact enough to generate consistent replica 
pressure fields. 

B. WEAKNESS OF THE SIMULATION 

In order to deal with reasonable computer times, only 
the cases of one or two unknown parameters were studied. 
For the same reason, the actual acoustic pressure field was 
compared with only a few replica fields. This limitation 
leads to moderate resolution in the results which could 
easily be improved by the generation of more replica fields. 

The weakness in modeling the noise field has already 
been emphasized in Chapter III. A. Further simulations 
should be done with a noise correlation matrix, KN, that has 
actually been obtained from measurements at sea. Further 
work using a more realistic parameterization of the ocean 
should be done, e.g., with empirical orthogonal functions 
(EOF) and using the technique to estimate EOF coefficients. 
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APPENDIX 



FORTRAN 77 PROGRAM USED IN THE SIMULATION 



PROGRAM ESTIMA 

ACOUSTIC TOMOGRAPHY USING MATCHED FIELD PROCESSING 

THIS PROGRAM DRAWS THREE KINDS OF DETECTION FACTOR IN 3D 
IT COMPARES BUCKER , BARTLETT AND MLM METHODS 

THE NOISF. IS I NT RODUCED BY ITS CORRELATION MATRIX 

THE NOISE MATRIX IS FROTORTIONAL TO IDENTITY MATRIX IN CASE OF 
UNCORRELATED NOISE BUT IS IN GENERAL HERMIT! AN MATRIX 

THE EXACT (OR MEASURED) PHASES AND MAGNITUDES ARE READ IN THE 
FILE EXACT DATA 

THE COMPUTED PARAMETERS ARE READ ON THE FILE NEAR DATA 

THE MODEL ALLOWS A MAXIMUM OF 20 RECEIVERS DUE TO THE USE OF THE 
PARABOLIC EQUATION FROM THE EARG PACKAGE 

KS IS THE COMPLEX MATRIX OF MEASURED PARAMETERS 
A IS THE MATRIX OF GUESSED PARAMETERS WE WANT TO MATCH 
NR IS THE NUMBER OF RECEIVERS (MAXI MUM 20) 

NF IS THE COMMON NUMBER OF POSSIBLE VALUES FOR 2 UNKNOWNS 
THE NUMBER OF FIELDS WE MATCH IS ACTUALLY NF*NF 



REAL rPIIASE (20) , PMAGNl (20).X(li),Y(ll ) . DF( 1 1 , 1 1 ) . DFCONT( 1 1 , I 1) 
REAL BART (11,11), FMLM (11.1 1 ) . MCONT( 1 1 . 1 1 ) . BCONT( I 1,11) 

COMPLEX A (20, 20) .KS(20,20) ,CDF( 1 I , 1 1 ) , KTT20 , 20) 

COMPLEX W(20 KCSUM, DETERM 
REAL KN(20 , 20) , NORM 

99 DEFINES THE FILE OF MEASURED ( EXACT) DATA 
98 DEFINES THE FILE OF DATA WE WANT TO MATCH 
NR=20 
NF=1 1 

READ THE MEASURED VALUES IN FILE EXACT DATA 
DO 1 1=1. NR 

READ (99 , 600) PMAGNl (I ) , PrilASE ( I ) 

CONTINUE 

FORMAT ( El 1 . 4 , 2X ,F7 . 4) 



999 



COMPUTE THE ELEMENTS OF MATRIX KS 
DO 2 J=-1.NR 
DO 3 K=1 ,NR 

KS(J,K)=PMAGNI (J) ,V PMAGNI (K) ,v EXP(CMPLX(0 . ,PPIIASE(J)- 
k PPIIASE (K) ) ) 

CONTINUE 

CONTINUE 
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Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr ->V Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr *>V Vr Vr Vr Vc Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr V r Vr Vr Vr Vr Vr V r Vr Vr V c Vr Vr Vr V r 



C DEFINE THE CORRELATION MATRIX KN OF THE CORRELATED NOISE 
PRINT*. VALUE OF SIGMA2 7 
READ* SIGMA2 

PRINT*. VALUE OF BETA ?' 

READ*, BETA 
DO 4 J=],NR 
DO 5 1=1. NR 

KN(I ,J)=SIGMA2*EXP(-BETA*ABS(I-J)) 

5 CONTINUE 

4 CONTINUE 



Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr >Y Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr 



c 

C COMPUTE THE TOTAL COVARIANCE MATRIX INCLUDING THE NOISE 
C BY ADDING THE NOISE AND THE SIGNAL CORRELATION MATRICES 
DO 6 J=1.NR 
DO 7 1=1, NR 
KT ( I , J ) =KS ( I , J ) +KN ( I ,J) 

7 CONTINUE 

6 CONTINUE 

C 
C 

C BEGINNING OF MAIN LOOP FOR THE NF RUNS 



C 



8 

C 

C 



10 

9 



DO 101 INDEX=],NF 
DO 100 JNI)EX=1,NF 

READ PARAMETERS OF EXPECTED FIRM) IN FILE NEAR DATA 
DO 8 1=1, NR 

READ(98 .600) PMAGNI ( I ) , PPIIASE ( l ) 

CONTINUE 

COMPUTE ELEMENTS OF MATRIX A 
DO 9 J=1,NR 
DO 10 K=1,NR 

A(J , K)=FMAGNI (J)*PMAGNI (K)*EXP (CMPLX (0 . , PPHASE(J)- 
FPIIASF, (K) ) ) 

CONTINUE 

CONTINUE 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C BUCKER DETECTION FACTOR C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 



12 

11 

C 



CDF ( INDEX , JNDEX )=CMPLX(0 . , 0 . ) 

DO 11 J=l , NR 
DO 12 K= 1 , NR 

I F(K . NE . J ) CDF( INDEX, JNDEX )- 
v A ( J . K)*CONJG (KT( J , K) ) 

CONTINUE 
CONTINUE 



‘CDF ( INDEX , JNDEX )+ 



C CI)F(- - , IS ACTUALLY REAL DUE TO PROPERTY OF MATRICES 

C A AND KT HENCE WE KEEP ONLY THE REAL PART OF IT 
DF( INDEX, JNDEX)=REAL(CDF( INDEX, JNDEX)) 

C 

C NORMALIZE THE DF FACTOR 
FACT0R=0 . 



FACT=0 . 



14 

13 



DO 13 J=1,NR 
DO 14 K=1 ,NR 

I F (K . NE . J ) FACTOR=FACTOR 1A ( J . K) *CONJG ( A ( J , K ) ) 
I F ( K . NE . J ) FACT=FACT+ KT ( J , K ) *CON JG ( KT (J , K) ) 
CONTINUE 
CONTINUE 



86 



FACTOR=SQRT( FACTOR) 

FACT-SQR i'(FACT) 

DF ( INDEX ,JNDEX)=DF (INDEX , JNUF.X )/ FACT/ FACTOR 
l)F (INDEX , JNDEX)=ARS (DF( INDEX , jflDEX) ) 

C 

100 CONTINUE 

101 CONTINUE 
C 

C NORMALIZE 13 Y THE LARGEST VALUE 

BHAX=0. 

00 67 1=1, NF 
DO 68 J=1 , NF 

1 F ( I)F ( 1 , J) .GE.BMAX) BMAX=DF( t , J) 

68 CONTINUE 

67 CONTINUE 

C 

DO 69 1=1, NF 
DO 70 J= 1 , NF 

DF( J , J )=DF ( I , J )/l3NAX 
70 CONTINUE 

69 CONTINUE 



C 

C DISSPLAY THE MATRIX OF NORMALIZED DETECTION FACTOR 

PRINT* , BUCKER FACTOR MATRIX 
DO 15 1=1, NF 

WRITE (6 ,777)(DF(I,J),J=1 , NF ) 

777 FORMAT ( 1 1 ( 2X , FA . 2 ) ) 

15 CONTINUE 

C 

cccccccccccccccccccccccccc 

C BARTLETT ESTIMATOR C 

CCCCCCCCCCCCCCCCCCCCCCCCCC 
REWIND 98 



C 



C 

c 



50 

C 

C 



16 

C 



17 

C 

C 



19 

18 



BEGINNING OF MAIN LOOP FOR IIIE NF RUNS 
DO 201 INDEX= I , NF 
DO 200 JNDEX=I , NF 

READ PARAMETERS OF REPLICA FIELD 
DO 50 1=1, NR 

READ (98 ,600) TMAGN I ( I ) , PPIIASE ( 1 ) 

CONTINUE 

DETERMINE THE NORMALIZED COMPLEX VECTOR W 
N0RM=0 . 

DO 16 1=1, NR 

NORM=NORM+PMAGNI (I)*+2 
CONTINUE 
NORM=SQRT(NORM) 

DO 17 1=1. NR 

W ( 1 ) = ( 1 . /NORM) *FMAGNI (I)* EXP(CMPLX(0. , PFIIASE( I ) ) ) 
CONTINUE 

COMPUTE THE BARTLEIT FACTOR 
CSUM=(0. ,0. ) 

DO 18 1=1, NR 
DO 19 J=1,NR 

CSUM=CSUM+CONJG (W ( I ) )*KT( I , J )*W ( J ) 

CONTINUE 

CONTINUE 

BART ( INDEX , JNDEX)=REAL(CSUM) 
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c 

200 

201 

C 



C 



21 

20 

C 



23 

22 

C 

C 



24 



CONTINUE 

CONTINUE 

NORMALIZATION BY THE HIGHEST VALUE 
BARMAX=0 . 

DO 20 1=1, NF 

DO 21 J=1,NF 

IF(BART(I,J) .GE. BARMAX) DARMAX=DART(I ,J) 
CONTINUE 
CONTINUE 

DO 22 1=1 .NF 

DO 23 J=1 ,NF 

BART( I , J)=BART(I ,J)/BARNAX 
CONTINUE 
CONTINUE 

DISSPLAY THE MATRIX OF NORMALIZED DETECTION FACTOR 
PRINT*. BARTLETT FACTORS MATRIX 
DO 24 1=1, NF 

WRITE (6, 7 7 7 ) (BAR 1(1,. I) , J=1 ,NF) 

CONTINUE. 



C 

CCCCCCCCCCCCCCCCCCCCC 
C MLM ESTIMATOR C 
CCCCCCCCCCCCCCCCCCCCC 



REWIND 98 



C 

C INVERT THE CORRELATION MATRIX KT 
CALL CM I RI N(NR , KT , NR , DETERM) 

C 

C BEGINNING OF MAIN LOOP FOR THE NF RUNS 
DO 301 JNDEX=1 ,NF 
DO 300 JNDEX=1 , NF 
C 

C READ PARAMETERS OF REPLICA FIELD 
DO 51 1=1, NR 

REAI)(98 , 600) PMAGNI ( I ) , PPHASE ( I ) 

51 CONTINUE 

C 

C DETERMINE THE NORMALIZED COMPLEX VECTOR W 

NORM-O . 

DO j-i MR 



26 

C 



27 

C 

C 



29 

28 

C 

300 

301 



NORM=NORMf PMAGNI ( I )**2 
CONTINUE 
NORM=SQRT ( NORN ) 

DO 27 1=1 NR 

W (I )-( 1 ■ /NORM)*PMAGNI (I)* EXP(CMPLX(0. , PPHASE ( I ))) 
CONTINUE 

COMPUTE THE MLM FACTOR 
CSUM=(0. ,0. ) 

DO 28 1=1. NR 
DO 29 J=1,NR 

CSUM=CSUM+CONJG(W(I))*KT(I ,J)*W(J) 

CONTINUE 

CONTINUE 

FMLM(INDEX,JNDEX)=REAL(l./CSUM) 

CONTINUE 

CONTINUE 
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c 

c 



31 

30 

C 



33 

32 

C 

C 



34 



NORMALIZATION BY THE HIGHEST VALUE 
FMLMAX=0 . 

DO 30 1=1, NF 

DO 31 J=1 ,NF 

IF(FMLM(I , J) . GE . FNLMAX) FMLMAX=FMLM(I , J) 
CONTINUE 
CONTINUE 

DO 32 1=1, NF 

DO 33 J=1,NF 

FULM(I , J)=FMLM(I , J)/FMLUAX 
CONTINUE 
CONTINUE 

DISSFLAY THE MATRIX OF NORMALIZED DETECTION FACTOR 
PRINT* , MLM FACTORS MATRIX' 

DO 34 1=1, NF 

WRITE (6, 777) (FMLM ( 1 , J) ,J=1,NF) 

CONTINUE 



C 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C COMPUTE THE SPREADING FACTOR OF THE ESTIMATORS C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccc 

WDF=0 . 

WB ART-0 . 

WFMLM=0. 

DO 650 1=1, NF 
DO 651 J=1,NF 

WDF=WDF I DF ( I , J)* (44 . I ] -40 . )**2* (54 I J -60 )**2 
WBART-WBAR IT BART( I , .J)* (44 . 1 1 - 40 . )**2* (54tJ -60)**2 
WFMLM=WFMLMIFMLM(I , J)* (44 . T I -40 . )**2*(54 . KJ-60 . )**2 
651 CONTINUE 

650 CONTINUE 
C 
C 

PRINT*, WDF= ,VDF 
TRINT* , WBART= .WBAPT 
PRINT*,' WFMLM= ' ,WFMI.M 
C 
C 

END 
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